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Abstract 

In this paper we extend the classical method of lattice dynamics to defective crystals with partial symmetries. 
We start by a nominal defect configuration and first relax it statically. Having the static equilibrium configu- 
ration, we use a quasiharmonic lattice dynamics approach to approximate the free energy. Finally, the defect 
structure at a finite temperature is obtained by minimizing the approximate Helmholtz free energy. For 
higher temperatures we take the relaxed configuration at a lower temperature as the reference configuration. 
This method can be used to semi-analytically study the structure of defects at low but non-zero tempera- 
tures, where molecular dynamics cannot be used. As an example, wc obtain the finite temperature structure 
of two 180° domain walls in a 2-D lattice of interacting dipoles. We dynamically relax both the position 
and polarization vectors. In particular, we show that increasing temperature the domain wall thicknesses 
increase. 
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1. Introduction 

Although it has been recognized that defects play an important role in nanostructurcd materials, the 
fundamental understanding of how defects alter the material properties is not satisfactory. The link between 
defects and the macroscopic behavior of materials remains a challenging problem. Classical mechanics of 
defects that studies materials with microscale defects is based on continuum theories with phenomcnolog- 
ical constitutive relations. In the nanoscale, the continuum quantities such as stress and strain become ill 
defined. In addition, due to size effects, to study defects in nano-structured materials, non-classical solu- 
tions of defect fields is necessary [^^j- The application of continuum mechanics to small-scale problems is 



problematic; atomistic numerical methods such as ab initio calculations [43|,|48[, Molecular Dynamics (MD) 



simulations [21|, l25| and Monte Carlo (MC) simulations [45|, |70| can be used for nanoscale mechanical anal- 
yses. However, the application of these methods is largely restricted by the size limit and the periodicity 
requirements. Current ab initio techniques are unable of handling systems with more than a few hundred 
atoms. Molecular dynamics simulations can model larger systems, however, MD is based on equations of 
classical mechanics and thus cannot be used for low temperatures, where quantum effects are dominant. 
Engineering with very small structures requires the ability to solve inverse problems and this cannot be 
achieved through purely numerical methods. What is ideally needed is a systematic method of analysis of 
solids with defects that is capable of treating finite temperature effects. 

The only analytic/semi-analytic method for solving zero-temperature defect problems in the lattice scale 
is the method of lattice statics. The method of lattice statics was introduced in [l^, E|. This method has 
been used for point defects [13, [3, for cracks lol [ill. [23 . and also for dislocations [J 11, 12, 4^ 55, 6S 



More details and history can be found in [34S [H JSiHOjlii , 4^ 5^ 63| and references therein. Lattice statics 
is based on energy minimization and cannot be used at finite temperatures. The other restriction of most 
lattice statics calculations is the harmonic approximation, which can be too crude close to defects. Recently, 
motivated by applications in ferroelectrics, we developed a general theory of anharmonic lattice statics 
capable of semi-analytic modeling of different defective crystals governed by different types of interatomic 
potentials [28[ . [68[ . [g^. At finite temperatures, the use of quantum mechanics-based lattice dynamics is 
necessary. Unfortunately, lattice dynamics has mostly been used for perfect crystals and for understanding 
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their thermodynamic properties [i,i,|3i,|3i,|4i,[5l|,|63. There is not much in the hterature on corrections 
for anharmonic effects and systematic solution techniques for defective crystals. Some of these issues will be 
addressed in this paper. 

In order to accurately predict the mechanical properties of nanosize devices one would need to take into 
account the effect of finite temperatures. It should be mentioned that most multiscale methods so far have 
been formulated for T — calculations. An example is the quasi-continuum method [49LJ53. However, 
recently there have been several attempts in extending this method for finite temperaturesjJQ IS 58 1. As 
Forsblom, et al. [l9j mention, very little is known about the vibrational properties of defects in crystalline 
solids. Sanati and Esetreicher 5J] showed the importance of vibrational effects in semi-conductors and the 
necessity of free energy calculations. Lattice dynamics 0, [Flj has been ignored with the exception of some 
very recent works 160|. As examples of finite-temperature defect solutions we can mention Taylor, et al. 
[59I 61 1 w ho discuss quasiharmonic lattice dynamics for three-body interactions in bulk crystals. Taylor, 
et al. [62| consider a slab, i.e., a system that is periodic only in two directions. They basically consider 
a supercell that is repeated in the plane periodically. As Allan, et al. P, Q conclude, a combination of 
quasiharmonic lattice dynamics, molecular dynamics, Monte Carlo simulations and ab initio calculations 
should be used in real applications. However, at this time there is no systematic method of lattice dynamics 
for thermodynamic analysis of defective systems that is also capable of capturing the anharmonic effects. We 
should mention that in many materials systems lattice dynamics is a valid approximation up to two-third 
of the bulk melting temperature but it turns out that harmonic approximation may not be adequate for 
free energy calculations of defects at high temperatures (see 18 1 for discussions on Cu). Hansen, et al. ^ 
show that for Al surfaces above the Debye temperature quasiharmonic lattice dynamic approximation starts 



to fail. Zhao, et al. [Tlj show that quasiharmonic lattice dynamics accurately predicts the thermodynamic 



properties of silicon for temperatures up to 800 K. In this paper, we are interested in low temperatures 
where MD fails while quasi-harmonic lattice dynamics is a good approximation. 

For understanding defect structures the main quantity of interest is the Helmholtz free energy. Free 
energy is an important thermodynamic function that determines the relative phase stability and can be used 
to generate other thermodynamic functions. In quasiharmonic lattice dynamics, for a system of n atoms, 
free energy is computed by diagonalizing a Sri x 3n matrix that is obtained by quadratizing the Hamiltonian 
about a given static equilibrium configuration. Using similar ideas, for a perfect crystal with a unit cell with 
TV atoms, one can compute the free energy by diagonalizing a 3A^ x 3A'^ matrix in the reciprocal space. In 
the local quasiharmonic approximation one assumes that atoms vibrate independently and thus all is needed 
for calculation of free energy is to diagonalize n 3 x 3 matrices [38j (see Rickman and LeSar 53 1 for a recent 
review of the existing methods for free energy calculations). These will be discussed in more detail in §3. 

In this paper, we propose a theoretical framework of quasi-harmonic lattice dynamics to address the 
mechanics of defects in crystalline solids at low but finite temperatures. The main ideas are summarized 
as follows. We think of a defective lattice problem as a discrete deformation of a collection of atoms 
to a discrete current configuration. The lattice atoms are assumed to interact through some interatomic 
potentials. At finite temperatures, the equilibrium positions of the atoms are not the same as their static 
equilibrium (T = 0) positions; the lattice atoms undergo thermal vibrations. The potential and Helmholtz 
free energies of the lattice are taken as discrete functionals of the discrete deformation mapping. For finite 
temperature equilibrium problems, the discrete nonlinear governing equations are linearized about a reference 
configuration. The finite-temperature equilibrium configuration of the defective lattice can then be obtained 
semi-analytically. For finite temperature dynamic problems, the Euler-Lagrange equations of motion of the 
lattice are casted into a system of ordinary differential equations by superimposing the phonon modes. We 
should emphasize that our method of lattice dynamics is not restricted to finite systems; defects in infinite 
lattices can be analyzed semi-analytically. The only restriction is the use of interatomic potentials. 

This paper is structured as follows. In §2 we briefly review the theory of anharmonic lattice statics 
presented in [g^ and [6^ . We then present an overview of the basic ideas of the method of lattice dynamics 
for both finite and infinite atomic systems in §3. This follows by an extension of these ideas to defective 
crystals with partial symmetries. In §4 we formulate the lattice dynamics governing equations for a 2-D 
lattice of dipoles with both short and long-range interactions. In §5 we study the temperature dependence 
of the structure of two 180° domain walls in the dipole lattice. Conclusions are given in §6. 
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2. Anharmonic Lattice Statics 



Consider a collection of atoms £ with the current configuration {x'j^^^ C M". Assuming that there 
is a discrete field of body forces {F'}ig£, a necessary condition for the current position {x*}ig£ to be in 
static equilibrium is —■§§[ + = 0, \/ i € C, where £ is the total static energy and is a function of the 
atomic positions. These discrete governing equations are highly nonlinear. In order to obtain semi-analytical 
solutions, we first linearize the governing equations with respect to a reference configuration Bq = {xg}ig£ 
[Hsj . We leave the reference configuration unspecified; at this point it would be enough to know that we 



usually choose the reference configuration to be a nominal defect configuration [28|, |68|, |69[ . 

Taylor expansion of the governing equations for an atom i about the reference configuration Bq = {xQ}ig£ 
reads 

^(^o)-^(^o).(x-x^)-Ea^(^o).(x^--x^„)-... + F'^^0. (1) 

Ignoring terms that arc quadratic and higher in {x^ — Xg}, we obtain 

^(^o).(x-x^) + E^(^o).(x-x^„) = -^(Z^„) + r y.eC. (2) 

jec 

Here, { — i^o)}i^c the discrete field of unbalanced forces. 

Defective Crystals and Symmetry Reduction. . In many defective crystals one can simplify the calculations by 
exploiting symmetries. A defect, by definition, is anything that breaks the translation invariance symmetry 
of the crystal. However, it may happen that a given defect does not affect the translation invariance of 
the crystal in one or two directions. With this idea, one can classify defective crystals into three groups: 
(i) with 1-D symmetry reduction, (ii) with 2-D symmetry reduction and (iii) with no symmetry reduction. 



Examples of (i), (ii) and (iii) are free surfaces, dislocations, and point defects, respectively j68l |. Assume 
that the defective crystal £ has a 1-D symmetry reduction, i.e. it can be partitioned into two-dimensional 
equivalence classes as follows 

N 



£ = U U (3) 



a6Z/=l 



where Si a is the equivalence class of all the atoms of type / and index a (see [68| and [28| for more details) . 
Here, we assume that £ is a multilattice of A'^ simple lattices. For a free surface, for example, each equivalence 
class is a set of atoms lying on a plane parallel to the free surface. Using this partitioning for i ~ la one 
can write 

S ^ ■ - ^ E' E E ■ (-" - 4') . (4) 

where the prime on the first sum means that the term J (3 — la is omitted. The linearized discrete governing 
equations are then written as [Hsf 

N / N 

E' E K/aJ^U^^ + - E' E I U^" = f/a, (5) 

^ez j=i y ^ez J=i 

where 

^i-JP= E d^w;^^^^^^ f/« = -a;^(So) + F,„, u-^^=x-^^-Xg^^=x^-x;^ Vje5,;,. (6) 
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The governing equations in terms of unit cell displacement vector Uq — (u^, u^)^ can be written as 

^A^(a)U„+^ =F„ aeZ, (7) 



where Ap{a) G M.^'^'^^'^, Uq,,Fq £ M'^^. This is a linear vector-valued ordinary difference equation with 
variable coefficient matrices. The unit cell force vectors and the unit cell stiffness matrices are defined as 



/ Kiq1/3 Kia2/J 
K2QI/3 K2a2/3 



a,/3 e Z. 



(8) 



Note that, in general, need not be symmetric l68i . The resulting system of difference equations can be 
solved directly or using discrete Fourier transform |68| . 

Hessian Matrix for the Bulk Crystal. . A bulk crystal is a defective crystal with a 0-D symmetry reduction. 
Governing equations for atom / in the unit cell n = read — + F/ = 0, / = 1,...,A^. Linearization 
about Bq = {X-^} yields 



Note that 



d£ 



3^1 



^{Bo)+Fi /=l,...,iV. (9) 



E ^(^0) . (x^- - X^) ^ ± E 7^(^o) . (x^ - X^) + E -^(^0) . (x^- - X^- 



j£C ,/=iie£j 
We also know that because of translation invariance of the potential 



jeCi 
3^1 



(10) 



Therefore, the linearized governing equations can be written as 



(11) 



N ( ^ \ 

J=l \ J=l / 



/ = l,...,iV, 



(12) 



where 



K 



1.1 



E 



(Bo) 



d£ 
ax^ 



The Hessian matrix of the bulk crystal is defined as 



H 



/ Kn K12 
K21 K22 



J ,J 



Kin \ 
K2JV 



X"' = x^ - X^' V J € Cj. (13) 



(14) 



where Kj/ = K/j. Stability of the bulk crystal dictates H to be positive-semidefinite with three zero 
eigenvalues. In the case of a defective crystal, one can look at a sequence of sublattices containing the defect 
and calculate the corresponding sequence of Hessians. 
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3. Method of Quasi-Harmonic Lattice Dynamics 



At a finite temperature T (constant volume) thermodynamic stability is governed by Helmholtz free 
energy F = E ~ TS. In principle, F is well-defined in the setting of statistical mechanics. Quantum- 
mcchanic ally calculated energy levels E{i) for different microscopic states can be used to obtain the partition 
function [33|, |66| 



Q = ^ exp 



(15) 



where kg is Boltzman's constant. Finally F — —kBTlnQ (see the appendix). However, one should note that 
the phase space is astronomically large even for a finite system. Usually, in practical problems, molecular 
dynamics and Monte Carlo simulations, coupled with thermodynamic integration techniques, reduce the 
complexity of the free energy calculations. For low to moderately high temperatures, quantum treatment of 
lattice vibrations in the harmonic approximation provides a reliable description of thermodynamic properties 
[4^ . In the following we review the classical formulation of lattice dynamics first for a finite collection of 
atoms and then for bulk crystals. 

3.1. Finite Systems 

For a finite system of atoms suppose B = {X^j^^^ is the static equilibrium configuration, i.e. 
0, y i & C. Hamiltonian of this collection is written as 



as I 



Now denoting the thermal displacements by u* = x* — X* potential energy of the system is written as 



Or 



£:(x)=£:(X) + V*u + o(|u|2), 



where $ is the matrix of force constants. The Hamiltonian is approximated by 

■H(x) ^ £{X) + iu^^u + iu"^Mu, 
where M is the diagonal mass matrix. Let us denote the matrix of eigenvectors of $ by [/. and write 

n{^) = £(X) + iqTAq + iq^Mq, 



(16) 



(17) 



(18) 



(19) 



(20) 



where q = U^u is the vector of normal displacements and A = diag(Ai, X^n) is the diagonal matrix of 
eigenvalues of This is now a set of SN independent harmonic oscillators. Solving Schrondinger's equation 
gives the energy levels of the ith oscillator as 44 1 



En 



: £:r(X) + [n + -]hujr n = 0,l, r = 1, SA^, 



(21) 



where Wr = t<J,- ^{X*} .^^j = \/ A,, / . The free energy is then written as Q 

3JV 00 

^{{^'}.eC'T) = -fcsT^lnJ^exp 



r—1 n— 



knT 



3N 3JV 



1 — exp 



huir 



(22) 
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Here it should be noted that we have considered a time-independent Hamiltonian, which can be regarded 
as a first-order approximation for some problems. Assume that Hamiltonian iJ of a system contains a time- 
dependent parameter t{t), say a time-dependent external force. If the time variation of f (t) is slow and does 
not cause a large variation of in a time interval of the same order as the natural period of the system with 
constant f, then this approximation is valid 1471 1 . otherwise one should consider time-dependent harmonic 



391 |42l. In such 



then this approximation is valid 
oscillator systems. This can be the case for various quantum mechanical systems 
situations one should obtain the solution of Schrondinger's equation for a time-dependent forced harmonic 
oscillator and as a result, energy levels would depend on the forcing terms too. As an example, Meyer [42| 
investigated energy propagation in a one-dimensional finite lattice with a time-dependent driving forces by 
solving the corresponding forced Schrondinger's equation. We also mention that the above formula for the 
free energy is based on the quasiharmonic approximation. As temperature increases such an approximation 
may become invalid for some materials 37| and therefore one would need to consider anharmonic effects. 



To include anharmonic terms in the free energy relation, anharmonic perturbation theory can be used by 
choosing the quasiharmonic state as the unperturbed state and the perturbation is due to the terms higher 
than second order in the Taylor expansion of the potential energy (Hij . This way, one accounts for anharmonic 
coupling of the vibrational modes. 

As we discuss in the appendix, to obtain the optimum positions of atoms at a constant temperature T 
one should minimize the free energy with respect to all the geometrical variables {X*}^^^ 
the governing equations are 



33|,l6^. Thus, 



d£ 



h duj. 



3N 



r=l 



1 



duj,. 



r=l CXp 



1 



ax» 



= 0. 



(23) 



To compute the derivatives of the eigenvalues, we use the method developed by Kantorovich [27|. Consider 
the expansion of the elements of the dynamaical matrix $ = [^afi] about a configuration B: 



iec 



E^(^)-(x-xo 



iec 



a,/3= 1, 



,3N. 



(24) 



If the eigenvectors of $ are normalized to unity, the perturbation expansion of eigenvalues would be [27 

3JV 



(25) 



iec a, 13 = 1 



where * denotes conjugate transpose and U = [Uap] is the matrix of eigenvectors of $ = [^'a^g], which are 
normalized to unity. Since higher order terms in the above expansion contain (x* — X')" with n G N > 2, 
all of them vanish for calculating the first derivatives of eigenvalues at x' ~ X' . Hence, wc can write 



dXr 
9x* 



3N 



_ d\r_ _ ^ d^afS rT 



a, 13=1 



(9X* 



and therefore 



dujr 



1 



2nirUJr 



3W 

1^ ^o^r-g^Upr- 



a, 13=1 



(26) 



(27) 



For minimizing the free energy, depending on the chosen numerical method, one may need the second 
derivatives of the eigenvalues as well. We can extend the above procedure and consider higher order terms 
to obtain higher order derivatives. The numerical method used in this paper for minimizing the free energy 
will be discussed in detail in the sequel. 



3.2. Perfect Crystals 

Let us reformulate the classical theory of lattice dynamics [1,11, 44 1 in our notation for a perfect crystal. 
This will make the formulation for defective crystals clearer. Let us assume that we arc given a multi-lattice 
C with simple sublattices. 



i.e. C = Ci- Let us denote the equilibrium position of i e £ by X* 
d 



9x* 



(28) 
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Atoms of the multi-lattice move from this equilibrium configuration due to thermal vibrations. Let us denote 
the dynamic position of atom z S £ by x* = We now look for a wave-like solution of the following 

form for i £ Cj 

:= X* - X' = ^U^(k) eK'^-^'--^'')*), 



/ITLj 



(29) 



where i = \/— T, w(k) is the frequency at wave number k e B, B is the first Brillouin zone of the sublattices, 
and U'^ is the polarization vector. Note that we are assuming that to/ 7^ 00 Note also that the displacements 
x*(i) are time dependent and are deviations from the average temperature-dependent configuration X* = 
X*(T). 

Hamiltonian of this system has the following form 



tec 



2 ^ 

iec 



■i\2 



-£({x^} 



(30) 



Because of translation invariancc of energy, it would be enough to look at the equations of motion for the 
unit cell S Z^. These read to/x^ 



I^l,...,N. Note that 



TO/X = 



/TO/ 



■U^(k)c.(k)2 e'(k-x^- 



.(k)t 



(31) 



The idea of harmonic lattice dynamics is to linearize the forcing term, i.e., to look at the following linearized 
equations of motion. 



m/x^ = -y-^(s)uJ" = -y y -^(s)u^ i = i, 

jec J=i]eCj 



.N. 



Note that for j € Cj 



— U'^(k) e'('^-^' 



/TOJ 



Therefore, equations of motion read 



TV 



c^(k)2u^(k) = yD/,/(k)U'^(k) 



,7=1 



where 



D 



ij 



(32) 



(33) 



(34) 



(35) 



are the sub-dynamical matrices. The case I = J should be treated carefully. We know that as a result of 
translation invariance of energy 

jec 

Thus 



TO/ ^ — ' 9xJ9x^ 

jeCi 

3^1 



d'8 



rrij ^ — ' dx^dx^ 
3^1 



(37) 



Finally, the dynamical matrix of the bulk crystal is defined as 
D(k) 



/ Dn(k) Di2(k) ... DiAr(k) \ 

D2l(k) D22(k) ... D2Ar(k) 



V DAri(k) DjV2(k) 



B3Afx3Af 



(38) 



Dwjv(k) / 



^For shell potentials, for example, shells are massless and one obtains an effective dynamical matrix for cores as will be 
explained in the sequel. 
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Let us denote the 3N eigenvalues of D(k) by Ai(k), i = 1, 3N. It is a well-known fact that the dynamical 
matrix is Hermitian and hence all its eigenvalues are real. The crystal is stable if and only if > V i. 
Free energy of the unit cell is now written as 



3JV 3N r / h (IrW 



k 1=1 k i=l 



where LUi = v/AlEI and a finite sum over k-points is used to approximate the integral over the first Brillouin 
zone of the phonon density of states. The second term on the right-hand side is the zero-point energy and 
the last term is the vibrational entropy. For the optimum configuration {X-' }^.^^ at temperature T, we have 

dT d£ v-v^ f h (l 1 \aa;,2(k)l 



9XJ 

Here using the same procedure as in the pervious section, one can calculate the derivatives of the eigenvalues 
as follows 

= 2^ U^., (k) Up, (k) , (41) 

where U (k) = [C/a^ (k)] e RS^^s^ is the matrix of the eigenvectors of D(k) = [^"'^(k)], which are 
normalized to unity. 

3.3. Lattices with Massless Particles 

Let us next consider a lattice in which some particles are assumed to be massless. The best well-known 
model with this property is the so-called "shell model" Let us assume that the unit cell has N particles 
(ions), each composed of a core and a (massless) shell. The lattice C is partitioned as 

N 
7=1 

Position vectors of core and shell of ion i are denoted by x* and x*, respectively. Given a configuration 
{x*}^g^, equations of motion for the fundamental unit cell read 

^.x^.-^, 0^-^, /^1,...,A.. (43) 

Assuming that cores and shells arc at a static equilibrium configuration, equations of motion in the harmonic 
approximation read 

mii^i = _y y ^!^.ui_y y ^!^.ui, / = i,...,iv, (44) 

J=i jeC} c^ccx^ j^;^ oxsox^ 

N q2£ ^ Q2g 

^ ^ ^ ^ ^ pi 3 pi I ^ ^ ^ ^ '1 3 P) I ' ^ l,...,iV. ('^^) 

J=l jeC-;, CXcCX^ j^;^ jg£s OXsOXg 

Note that for j e Cj wc can write 

-^U'^(k) ei(k-x^e-(k)0^ = uf (k) e*('^-^'=-"('')0 k e B, (46) 



/mj 



^Note that this is consistent with Eq. I I21I I as wc arc using mass-reduced displacements. 
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where B is the first Brillouin zone of £J (or £f ). Tlius, and ean be simplified to read 



N 



N 



^ D- U,^(k) + ^ (k) = c.2(k)U^(k) / ^ 1, iV, 



N 



J=l 

AT 



^ DlSU'^(k) + J2 DljU,^(k) = 



j=i 



where 



D 



/ = l,...,iV, 



D 



I. J 



E 



,ik-(XJ,-X^ 



Eqs. (|T7|) and (|3S]) can be rewritten as 



DecU, + D,,U, = uj^Vc and U, = -D;/D,,U„ 



where 




Nl 




D 



IN 



D 



NN 



Finally, the effective dynamical problem for cores can be written as 

D(k)Ue(k) =w(k)2U,(k), 

where 



(47) 
(48) 



(49) 

(50) 

(51) 
(52) 
(53) 

(54) 
(55) 



D(k) - D,,(k) - D,,(k)D-i(k)D,,(k), 

is the effective dynamical matrix. Note that T)cs and T)sc are not Hermitian but DcsDj/Dsc is. 

The diagonal submatrices of D, i.e. Djj and Dff should be calculated considering the translation 
invariance of energy, namely 



Dec _ 
// — 



J_ \ ' ^ ^ ^ik-(xj-x5) ^ \ ^ 9 £ 

mi dxidxl mi ^ 5xJ(9x^' 



(56) 
(57) 



Denoting the 3A^ eigenvalues of D(k) by Ai(k) = wK'*)' energy of the unit cell is expressed as 



J-({X^„Xi},e£,T) =£({X^„Xi},e£) -nc.,(k) + fcBTln 



k 1=1 



1 — exp 



^i(k) 



(58) 
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Therefore, for the optimum configuration {X^jX^}^.^^ at temperature T we have 




dJ- d£ J " I ^ ^ 1 '-'^iK'^) [ /ar\\ 

where the derivatives of eigenvalues are given by 

Bum ^ y y. ^^^v,., (k) , MM ^ y (k) ^^^y,, (k) , (61) 

^ a,p — l ^ ^ a,p — l ^ 

where V (k) = [Vq^ (k)] e m3^><3^ is the matrix of the eigenvectors of D(k) = [Dap{}<i)], which are 
normalized to unity. 

3.4^. Defective Crystals 

Without loss of generality, let us consider a defective crystal with a 1-D symmetry reduction [HI], i.e. 

N 

^-^UU^JP- (62) 

J=l /36Z 

Note that j = J/3 means that the atom j is in the (5th equivalence class of the Jth sublattice. For this atom 
the thermal displacement vector is assumed to have the following form 



u-'" = -^U^^(k) ei(k-x^-"(k)*), k € B, 



/mj 

where B is the first Brillouin zone of Cj. Equations of motion in this case read 

N 



(63) 



a.(k)2u^-(k) = Y.Y. D/aJ/3(k)U^^(k), (64) 

j=i ,9ez 

where 

^....-^^^^^-•--"'^,im, (65) 

are the dynamical sub-matrices. The sub-matrices Tiiaia have the following simplified form 

D/./. = - E -''■'^"^^°'7r4S^(^)- (66) 

mi ^-^ ox'^oxJ 
i6£ 



Note that 



Thus 



- - E ^"■^"^""^"^7r4S^(^) - - E 7r4^(^)- (68) 

It is seen that for a defective crystal the dynamical matrix is infinite dimensional. 

As an approximation, similar to that presented in [s^ as the local quasiharmonic approximation, one 
can assume that given a unit cell, only a finite number of neighboring equivalence classes interact with its 
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thermal vibrations. One way of approximating the free energy would then be to consider vibrational effects 
in a finite region around the defect and study the convergence of the results as a function of the size of the 
finite region. For similar ideas see [1^, H^l, and 13 1. Here, we consider a finite number of equivalence classes, 
say ~C < a < C, around the defect and assume the temperature-dependent bulk configuration outside this 
region. As another approximation we assume that only a finite number of equivalence classes interact with 
a given equivalence class in calculating the dynamical matrix, i.e. we write 



m N 
a— — m 7—1 

where Ci is the neighboring set of atom i. Therefore, the linearized equations of motion read 

m N 

c.(k)2u^"(k) = ^D,<,j^(k)U^'^(k) a = -C,...,C. 



Defining 



l3 = -m .1=1 



u 



Nc 



we can write the equations of motion as follows 



(69) 



(70) 



(71) 



;(k)2u„(k)= A„(„+^)(k)U(„+^)(k), 



where 



Di 



Ql/3 



p3Afx3Af 



D 



Naif) 



D 



NaNp 



(72) 



(73) 



Now considering the finite classes around the defect, we can write the global equations of motion for the 
finite system as 

D(k)U(k) = a;(k)2u(k), (74) 

where 



U(k) 




, D(k) 



and 



>'(-C){-C) 
Aa^ 



-C)C 



"CC 



aMxM 



, M = 37Vx (2C+1), (75) 



(76) 



It is easy to show that AQ,^(k) = A^^(k), i.e. the dynamical matrix D(k) is Hcrmitian, and therefore has 
AI real eigenvalues. Note that the defective crystal is stable if and only if ujf > V i. 
Now we can write the free energy of the defective crystal as 



M , 



k i=l 



1 — cxp — 



Huji (k) 



(77) 



In the optimum configuration {X-'}^.^^ at a finite temperature T, we have 



ax. 



k i= 




(78) 
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where the derivatives of the eigenvalues are calculated as follows 

- 2^ U^r (k) gy^j Up, (k) , (79) 

a,/3 — 1 

where U (k) = \Uap (k)] G jj^/x a/ -g |-j^g matrix of the eigenvectors of D(k) = [Dap (k)], which are normalized 
to unity. 

3.5. Defect Structure at Finite Temperatures 

In the static case, given a configuration B'q = \ x'n > , one can calculate the energy and hence forces 

exactly, as the potential energy is calculated by some given empirical interatomic potentials. Suppose one 
starts with a reference configuration and solves for the following harmonic problem: 

^ {B',)-{^^^A)^-^{B',) VieC. (80) 



^ c'x'c'xJ ' ' ax 

This reference configuration could be some nominal (unrelaxed) configuration. Then one can modify the 
reference configuration and by modified Newton- Raphson iterations converge to an equilibrium configuration 
Bq = {xo}jg£ assuming that such a configuration exists [g^. In this configuration (Bq) = 0, \l i € C. 
Bo is now the starting configuration for lattice dynamics]! For a temperature T, the defective crystal is in 
thermal equilibrium if the free energy is minimized, i.e., if 

^(B)^0,y^eC. (81) 

Solving this problem one can modify the reference configuration and calculate the optimum configuration. 
This iteration would give a configuration that minimizes the harmonically calculated free energy. The next 
step then would be to correct for anharmonic effects in the vibrational frequencies. One way of doing this is 
to iteratively calculate the vibrational unbalanced forces using higher order terms in the Taylor expansion. 

There are many different optimization techniques to solve the unconstrained minimization problem ()8ip . 
Here we only consider two main methods that are usually more efficient, namely those that require only the 
gradient and those that require the gradient and the Hessian (sij . In problems in which the Hessian is avail- 
able, the Newton method is usually the most powerful. It is based on the following quadratic approximation 
near the current configuration 

^ (b'' + = J- (^''0 + VJ- (S^-) • + ^{d'V ■ H {B") -s' + o (l^V) , (82) 

" k ~k 

where 8 — B^^^ — B''. Now if wc differentiate the above formula with respect to ^ , we obtain Newton 
method for determining the next configuration B'^^^ ~ B'^ + S : 6'^ = — H^^ (S'^) • {B^^. Here in 
order to converge to a local minimum the Hessian must be positive definite. 

One can use a perturbation method to obtain the second derivatives of the free energy but as the dimen- 
sion of a defective crystal increases, calculation of these higher order derivatives may become numerically 
inefficient (g^l and so one may prefer to use those methods that do not require the second derivatives. One 
such method is the quasi-Newton method. The main idea behind this method is to start from a positive- 
definite approximation to the inverse Hessian and to modify this approximation in each iteration using the 
gradient vector of that step. Close to the local minimum, the approximate inverse Hessian approaches 
the true inverse Hessian and we would have the quadratic convergence of Newton method (H^ . There are 
different algorithms for generating the approximate inverse Hessian. One of the most well known is the 
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm js^ : 

A.+ .=A. + i;^-<^^^M^ + (AT.A-.A)u8u, (83) 



^If temperature is "large", one can start with equilibrium configuration of a lower temperature. This is what wc do in our 
numerical examples as will be discussed in the sequel. 



12 



where A* = (H* 



A = VJ-«+i - VF\ and 



u 



■ A 



A' • A 
AT-A'- A' 



(84) 



Calculating A*+^, one then should use A'+^ instead of H ^ to update the current configuration for the 

~ k 

next configuration B'^^^ ~ B'^ + 6 . If A*+^ is a poor approximation, then one n iay need to perform a 
linear search to refine i5'=+^ before starting the next iteration 52 1. As Taylor, et al. [60[ mention, since the 
dynamical contributions to the Hessian are usually small, one can use only the static part of the free energy 
£ to generate the first approximation to the Hessian of the free energy. Therefore, we propose the following 
quasiharmonic lattice dynamics algorithm based on the quasi-Newton method: 



Input data: Bq (or Bt-at for large T), T 
O Initialization 

t> = H.static\B=Bo 

[> Do until convergence is achieved 

> T>^ ^ TtiB^) 

> Calculate V J"'" 

[> Use A*-' to obtain B^^^ 

t> End Do 
O End 



4. Lattice Dynamic Analysis of a Defective Lattice of Point Dipoles 



In this section we consider a two-dimensional defective lattice of dipoles. Westhaus [67[ derived the 
normal mode frequencies for a 2-D rectangular lattice of point dipoles using the assumption that interacting 
dipoles have fixed length polarization vectors that can only rotate around fixed lattice sites. In this section, 
we relax these assumptions and in the next section will obtain the temperature-dependent structures of two 
180° domain walls. 

Consider a defective lattice of dipoles in which each lattice point represents a unit cell and the corre- 
sponding dipole is a measure of the distortion of the unit cell with respect to a high symmetry phase. Total 
energy of the lattice is assumed to have the following three parts 



£ ({x', P^I.e^) = ({x^ P'},e£) + f {{^'}^ec) + ({P'}.g£ 



(85) 



where, £'*'hort ga. .^^^^ ^^^^ dipole energy, short-range energy, and anisotropy energy, respectively. The 
dipole energy has the following form 



2 ^ 



P* • P^' 3P' • (x* - x^ ) P^' • (x* - 



E- 

^ 2a,; 



(86) 



where is the electric polarizability and is assumed to be a constant for each sublattice. For the sake of 
simplicity, we assume that polarizability is temperature independent. The short-range energy is modeled by 
a Lennard- Jones potential with the following form 



E 



short 



\ E 









V |x* -x^ l 




V|x» -xJ| / 



(87) 
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where for a multi-lattice with two sublattices Uij and eij take values in the sets {an, ai2, 022} and {en, ei2, £22}, 
respectively. The anisotropy energy quantifies the tendency of the lattice to remain in some energy wells 
and is assumed to have the following form 



= ^X^|P' -Pi|2 |P^-P 



2 

21 ■ 



This means that the dipoles prefer to have values in the set {Pi, P2}- 

Let S ~ ({X*,P*}ig£) be the equilibrium configuration (a local minimum of the energy), i.e. 

= yieC. (89) 



9P» 



It was shown in [68| how to find a static equilibrium equation starting from a reference configuration. We 
assume that this configuration is given and denote it by yB = {X*, P'}ig£. At a finite temperature T, ignoring 
the dipole inertia, Hamiltonian of this system can be written as 



H({x',P^},e£) ^ij^m^lxf +£:({x\P^},e£). (90) 

Equations of motion read 



Linearizing the equations of motion (j9ip about the equilibrium configuration, we obtain 
-m,x' = ^^'f ■ (B) (x' - X') + V ■ (B) (xJ - X^) 

where Si = C\ {i}. Note that 

^p-^ (B) = 2if^ (|P^ - Pip + |P* - P2n I + ^Ka (P' - Pi) ® (P^ - P2) 

+ AKa (P' - P2) ® (P' - Pi) + —I, (94) 

where I is the 2x2 identity matrix and (E) denotes tensor product. 

For a defective crystal with a 1-D symmetry reduction the set C can be partitioned as follows 

TV 

>C=|J|j£7a. (95) 

a£ZI=l 

Let us define u* = x' — X*, q* = P' — P*. Periodicity of the lattice allows us to write for i e Cia 

= ^Lu^"(k) e'(k-x'-<^(k)t)^ q» ^ Q^"(k) e'C'-^'-^C^)*)^ k e B. (96) 



/mj 



14 



Thus, Eq. for i ^ la can be simplified to read 

d'S , w 1 d'S 



.7=1 ^iez ie£j3 V 

V J=l ^gZ je£.;^ V 

where a prime on summations means that the term corresponding to Jf3 ~ la is excluded. Eq. can be 
rewritten as 

N N 
J=l I3£Z J=1/3£Z 

where 



Similarly, Eq. (|93p can be simplified to read 



Or 



where 



.1=1 pel, jeCjn V 

+ ^^.^ (^) Q^" W + E E E' ap?^ (^) e-(-^--^")Q^^(k) = 0. (100) 

N N 

E E D?:,7,(k)u^^(k) + E E Df.,,(k)Q'^^(k) = 0, (101) 

j=i /3ez J=i /3ez 



Df„,,(k) =^.,^,.^p^(S)+ E' ap?^(^)^"-^"^-"^°^- (102) 



We know that 166 



And 



Before proceeding any further, let us first look at dynamical matrix of the bulk lattice. 
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Dynamical Matrix for the Bulk Lattice.. In the case of the bulk lattice we have 

N 



c=\Jc,. 



1=1 



Periodicity of the lattice allows us to write for i G £/ 

= -Lu^(k) e'^k-^'-^C')*), = Q^(k) eHk-x'--(k)t)^ ^ g. 



/rni 

Thus, Eq. (|9T|) for i = / is simplified to read 

1 d^£ 



uikyv' (k) 



TO/ dx^dx. 

1 d^S 



, iB) U^(k) + E E' 77^7^ (^) e-^--^)U^(k) 



J=i ie£,7 

AT 



y/rnjrnj dx^dx^ 
J 1 



-(S)Q^(k) + E E -^7^^(S)e*''(^^-^^)Q^(k). 



J=i ie£j 



This can be rewritten as 



AT 



N 



;(k)2u^(k) = E D?j(k)U'^(k) + E D^5(k)Q'^(k), 



j=i 



j=i 



where 



Df5(k) = 5,A,^(fi) + E' ' 



mi dx^ dx^ 
1 



y/rnpmj dx^dx 



j{B)e 



ik-(XJ-X^) 



/TO/ ap^9x 

Similarly, Eq. ([9T|) is simplified to read 

1 d^s 



/m/ dV^dx 



N 



/mi dx^dV^ 



(S)U^(k) + E E' 77^7^1^ (^)e-(x-x^)u./(k) 



je£j 



(^)Q'(k) + E E'a^^^)^ 



J=i je£j 



ik-(X'-X^)Q^(k) = 0. 



Or 



where 



AT 



N 



E D?^(k)U'^(k) + E Df,(k)Q'^(k) = 0, 



j=i 



Dr/(k)=<5/,/- 



/TO/ dxWP^ 



fmj dx^dP^ 



Defining 











-1 
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the linearized equations of motion read 

D,,(k)U(k) + D,p(k)Q(k) = c^(k)2u(k), Dp,(k)U(k) + Dpp(k)Q(k) = 0, 

where 



(114) 





/ Dff . 


J\XX \ 




/ D^f . 


'-'in 
















I ^Wi ■ 


■ '-'nn / 






• ■ '-'nn 




f D?^ . 


'-'in ' 




f D?f . 


'-'in 


"^px 














\ ^Nl ■ 








■ '-'nn 



(115) 

1-)PP T^PP j 

Finally, the effective dynamical problem can be written as 

D(k)U(k) = w(k)2u(k), (116) 

where 

D(k) = D,,(k) - D,p(k)D;pi(k)Dp,(k), (117) 

is the effective dynamical matrix. Note that D(k) is Hermitian. Denoting the 2N eigenvalues of D(k) by 
Ai(k) = a;^(k), i — 1, 2N, free energy of the unit ceU is expressed as 



2W , 

^ ({X^ P^"},e£, T) = ({X^ P^I.e^) + ^ ^ j -fiw,(k) + fcsTln 



1 — exp 



huJi{k) 
knT 



(118) 



Therefore, for the optimum configuration |X^ , P^ }_^.^^ at temperature T we should have 



dJ- d£ 



2N 



8X3 0X3 



dT dS 



EE|2..(k) ^2 ■ exp(^)-l 



2N 



^l2-.(k)\2 expf^)_i 




= 0, 



(119) 
(120) 



where the derivatives of eigenvalues are given by 



9^1 (k) 

Q,/3 — 1 ct.p — 1 

where V (k) = [V^^ (k)] e ^2^x2^ ^j^^, matrix of the eigenvectors of D(k) = [^'^^(k)], with D^^n 
normalized to unity. 

Dynamical Matrix for the Defective Lattice. In the case of a defective lattice we consider interactions of 
order m, i.e., we write 

m N 
a — — rn I—l 

where £i is the neighboring set of the atom i. The equations of motion (|M1) and (|10ip become 

N m N m 

a.(k)2u^-(k) = E E D?S^^(k)U^''(k) + ^ ^Tcjpim'^i^), (123) 

.J=l f} = ~m .J=l l3 = -m 

N m N m 

= E E D?:,,(k)U'^^(k) + 5: J2 Df„.,(k)Q'^^(k). (124) 
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Defining 



r,2N 



u 



Na 



i2N 



we can write tlie equations of motion as follows 

m m 

a;(k)2u„(k) = ^ A-^„+^)(k)U(„+^)(k)+ ^ A^^^^^)(k)Q(„+^)(k), 



/3=-m 



/3 = -m 







E A:^o+^)(k)U(„+,)(k)+ ^ A^4^,)(k)Q(„+,)(k), 

/3 — — m (5 ——m 



(125) 

(126) 
(127) 



where 



A** — 



^NaNp 



eM2^><2Ar ^^^^^^^^ 



(128) 



Let us consider only a finite number of equivalence classes around the defect, i.e., we assume that — C < a < 
C . Therefore, the approximating finite system has the following governing equations 



D,,(k)U(k) + D,p(k)Q(k) = w(k)2u(k), 
Dp,(k)U(k)+Dpp(k)Q(k) =0, 



where 



U 



c 



D,4k) 



V 



U(k) = 



eM^^ Q(k) = 



Q c 
Qc 



"(-0)0 



A- |a-/3|<m. 



(129) 
(130) 



(131) 



, (132) 



where M — 2N x (2C + 1) and *,-k = x,p. Now the effective dynamical problem can be written as 

D(k)U(k) = a;(k)2U(k), (133) 

where 

D(k) = D,,(k) - D,p(k)D;pi(k)Dp,(k), (134) 
is the effective dynamical matrix. Note that D(k) is Hermitian and has M real eigenvalues. The free energy 



of the unit cell is expressed as 



k i=l 



1 — exp — 



^^i(k) 
knT 



(135) 



For the optimum structure {X-' , P-' }^.^^ at temperature T we have 



d£ 



M 



eeI 



dX^ ^ ^ I 2..,(k) \ 2 _ 1 



dT _ dE_ 



eeI 



k i=l 



2c^,(k) 2 



exp 




0, 



(136) 
(137) 
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where the derivatives of eigenvalues are given by 



a,/3— 1 a,^— 1 

where V (k) = [Vq^ (k)] G M^^^'*^ is the matrix of the eigenvectors of D(k) = [D^p (k)], with normalized 
to unity. 



5. Temperature-Dependent Structure of 180° Domain Walls in a 2-D Lattice of Dipoles 

To demonstrate the capabilities of our lattice dynamics technique, here we consider a simple example of 
180° domain walls shown in Fig. [T] In these 180° domain walls, polarization vector changes from — Po on 
the left side of the domain wall to Pq on the right side of the domain wall. We consider two types of domain 
walls: Type I and Type II. In Type I (the left configuration) the domain wall is not a crystallographic line, 
but it passes through some atoms in Type II (the right configuration). We are interested in the structure of 
the defective lattice close to the domain wall at a finite temperature T. In these examples, each equivalent 
class is a set of atoms lying on a line parallel to the domain wall, i.e., we have a defective crystal with a 
1-D symmetry reduction. The static configurations for Type I domain wall, Bq, was computed in jisj . Here 
we consider the static equilibrium configurations as the initial reference configurations. For index n £ Z 
in the reduced lattice (see Fig. [1]), the vectors of unknowns are U„,Q„ e R'^. Because of symmetry, we 
only consider the right half of the lattices and because the effective potential is highly localized [68[, for 
calculation of the stiffness matrices, we assume that a given unit cell interacts only with its nearest neighbor 
equivalence classes, i.e., we consider interactions of order m = 1. Note that this choice of m only affects 
the harmonic solutions; the final anharmonic solutions are not affected by this choice. For our numerical 
calculations we choose = 280 atoms in each equivalence class as the results are independent of N for 
larger N . Note that for force calculations we consider all the atoms within a specific cut-off radius Re- Here, 
we use Rc = 140a, where a is the lattice parameter in the nominal configuration. 



Symmetry 



Symmetry 



-1 1 



-1 



Figure 1: Reference configurations for the 180° domain walls in the 2-D lattice of dipoles, their symmetry reduction and their 
reduced lattices. Left Panel: Type I, Right Panel: Type II. 



For minimizing the free energy, first one should calculate the effective dynamical matrix according to Eq. 
(|134p . The calculations of this matrix for the two configurations are similar. For example, in configuration 
I due to symmetry we have U_i = —Up. Also we consider the temperature-dependent bulk configuration 
as the far- field condition, i.e., we assume = Uc for a > C + 1. Our numerical experiments show that 
choosing C = 35 would be enough to capture the structure of the atomic displacements near the defect, so 
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we use C = 35 in what follows. For the right half of the defective lattice we have 



where 5 = 2 (C 
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1), 

E^ 



D 



00 



D 



0(-l) 



and 



c 



D 



cc 



D 



C(C+1) 



pSxS 



(139) 



a;, p. 



(140) 



Now one can use the above matrices to calculate the effective dynamical matrix. Note that as a consequence 
of considering interaction of order m, the dynamical matrix will be sparse, i.e., only a small number of 
elements are nonzero. As the dimension of the system increases, sparsity can be very helpful in the numerical 
computations [52 1. 




Figure 2: Position and polarization displacements for Type I domain wall (T = 5) obtained by choosing different number of 
k-points (r) in the integration over the first Brillouin zone. 

As was mentioned earlier, we will consider only the static part of the free energy to build the Hessian for 
the initial iteration and then update the Hessian using the BFGS algorithm in each step. To calculate the 
gradient of the free energy we need the third derivatives of the potential energy. These can be calculated 
using following relation 



Q's ~ g-s Q's pp '-'xp'-'pp '-'pp '-'px '-'xp'-'pp Q„ ^ — ^ ,r . v-^^-^; 



To obtain these third derivatives one can use the translation invariance relations (|103p and (|104p to simplify 
the calculations. For example, we can write 

jec, 

where a prime means that we exclude j = i from the summation. 

The dimensionalized temperature T and dimensionalized mass fh correspond to the choice h = = 
and work with normalized . To obtain the static equilibrium configuration and also in dynamic 



* We select these values to be able to work with temperatures that are comparable with real temperature values. 
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calculations we use a = 1.0, Po = 1.0, e = 0.125 , Ka = 2.0 and to = lO''. In what follows convergence 
tolerance for VVJ^ • VJ^^ is 10^^. Using this value for convergence tolerance, solutions converge after ten 
to twenty iterations. In Fig. [5] we plot and qy for Type I domain wall and T = 5 for different number 
of k-points (r) in the first Brillouin zone. Here is the diplacement of the lattice with respect to the 
nominal configuration at temperature For numerical integrations over the first Brillouin zone we use the 
special points introduced in l46[. For the case r = 1 we set k = 0, i.e., we assume that all of the atoms in 
a particular equivalence class vibrate with the same phase. As can be seen in these figures, displacements 
converge quickly by selecting r = 7 k-points in the first Brillouin zone, so in what follows we set r = 7. 
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Figure 3: Position and polarization displacements of Type I domain wall with respect to the temperature-dependent nominal 
configurations. 
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Figure 4: Position and polarization displacements of Type II domain wall with respect to the temperature-dependent nominal 
configurations. 



Figs. [3] and 2] show the variations of displacements with temperature for the two domain walls. As 
temperature increases we cannot use the static equilibrium configuration as the reference configuration for 
calculating Hp. Instead, we use the equilibrium configuration of a smaller temperature to obtain Hq. Here, 
we use steps equal to AT ~ 5. In other words, for calculating the structure of a domain wall at T = 30, for 
example, we use the structure at T ~ 25 as the initial configuration. We see that the lattice statics solution 
and the lattice configuration at T = obtained by the free energy minimization have a small difference. 
Such differences are due to the zero-point motions; the lattice statics method ignores the quantum effects. 
It is a well known fact that zero-point motions can have significant effects in some systems 3l| . Note that 



'""Note that as temperature increases, lattice parameters change. A temperature-dependent nominal configuration is what is 
shown in Fig. [T]but with the bulk lattice parameters at that temperature. 
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polarization near the domain wall increases with temperature. Also as it is expected, the lattice expands by 
increasing the temperature. 

Only a few layers around the domain wall arc distorted; the rest of the lattice is displaced rigidly. As we 
see in Fig. [SJ the domain wall thickness for both configurations increases as temperature increases. In this 
figure Wt = wt/wq^ where wq is the domain wall thickness at T = 0. Note also that in this temperature 
range wt increases linearly with T. This qualitatively agrees with experimental observations for PbTiOs 
in the low temperature regime [l7|. Foeth, et al. 17[ observed that domain wall thickness increases with 



temperature. What they measured was an average domain wall thickness. Note that domain wall thickness 
cannot be defined uniquely very much like boundary layer thickness in fluid mechanics. Here, domain wall 
thickness is by definition the region that is affected by the domain wall, i.e. those layers that are distorted. 
One can use definitions like the 99%-thickness in fluid mechanics and define the domain wall thickness as 
the length of the region that has 99% of the far field rigid translation displacement. What is important is 
that no matter what definition is chosen, domain wall "thickness" increases by increasing temperature. 

Our calculations show that by increasing the mass of the atoms both position and polarization displace- 
ments decrease. However, variations of displacements with respect to mass is very small. For example, by 
increasing mass from ffi = lO"' to m = 10^ at T = 10, displacements decrease by less than 0.1%. 





Figure 5: Variation of the 180° domain wall thickness with temperature. 



6. Concluding Remarks 

In this paper we extended the classical method of lattice dynamics to defective crystals. The motivation 
for developing such a technique is to semi-analytically obtain the finite-temperature structure of defects in 
crystalline solids at low temperatures. Our technique exploits partial symmetries of defects. We worked out 
examples of defects in a 2-D lattice of interacting dipolcs. We obtained the finite-temperature structure of 
two 180° domain walls. We observed that using our simple model potential, increasing temperature domain 
walls thicken. This is in agreement with experimental results for ferroelectric domain walls in PbTiOa. This 
technique can be used for many physically important material systems. Extending the present calculations 
for 180° domain walls in PbTiOs will be the subject of a future work. 



Appendix A. The Ensemble Theories 

There are different ensemble theories for calculating the thermodynamical properties of systems from the 
statistical mechanics point of view. In this appendix, we consider micro canonical and canonical ensemble 
theories and discuss the relation between them. In particular, we will see that the free energy minimization 
discussed in this paper i s eq uivalent to finding the most probable energy at the given temperature. For more 
detailed discussions see |50l . 
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Appendix A.l. Micro Canonical Ensemble Theory 

Prom thcrniodynamical considerations, it is known that by specifying the hmited number of properties 
of a system, one can determine ah the other properties. In principle, any physical system, i.e., any macro 
system, consists of many smaller subsystems. Therefore, we can consider properties of each macro system 
as macrostates specified by the properties of these subsystems that are called microstates. Note that by 
a microstate we mean a set of values associated to each subsystem of a system. For example, consider an 
isolated system with energy E and volume V that consists of N non- interacting particles with energies Ci, 
i = 1,2, . . . ,N. Now each n-topple (ei, . . . , e^) satisfying 

N 

J2e, = E, (A.l) 

i=l 

would represent a microstate of this system. 

Obviously, there may exist several microstates that are associated to the same macrostate. Let ^{E, N, V) 
denote the number of microstates associated with the given macrostate {E, N, V). We assume that for an 
isolated system, (i) all microstate compatible with the given macrostates are equally probable, and (ii) 
equilibrium corresponds to the macrostate having the largest number of microstates. Let S and ks denote 
the entropy of a system and Boltzmann constant, respectively. Then one can show that the above two 
assumptions and setting 

S^kB\n{n), (A.2) 

yields the equality of temperatures for systems that are in thermodynamical equilibrium. Note that (|A.2[) 
provides the fundamental relation between thermodynamics and statistical mechanics. Once S is obtained, 
the derivation of other thermodynamical quantities would be a straight forward task. 

Appendix A.2. Canonical Ensemble Theory 

In practice, we never have an isolated system and even if we have such a system, it is hard to measure 
the total energy of the system. This means that it is more convenient to develop a statistical mechanics 
formalism that does not use E as an independent variable. It is relatively easy to control the temperature 
of a system, i.e. we can always put the system in contact with a heat bath at temperature T. Thus, it is 
natural to choose T instead of E. 

Let a system be in equilibrium with a heat bath at temperature T [^. In principle, the energy of the 
system at any instant of time can be equal to any energy level of the system. As a matter of fact, one can 
show that the probability of a system being in the energy level Pr is equal to 

^ gr CXp{-Er/kBT) ^ Qr CXp{-Er/kBT) 

E^9^^M-EJkBT) Q(r,T) ' ^ • ^ 

where we define the partition function of the system as 

Q(r,T) = ^5,exp(-S,/fcBr), (A.4) 

i 

and T denotes any other parameters that might govern the values of Er . Note that the summation goes over 
all energy levels of the system and gi denotes the degeneracy of the state Ei, i.e. the number of different 
states associated with the energy level Ei. Thus, one may write gi ~ Cl(Ei), where O comes from the previous 
formulation. Assuming the total energy of the system to be an average energy of the different states, i.e. 



E = J2PrEr, (A.5) 



one can show that the Helmholtz free energy T can be written as 

J-=-kBT\nQ. (A.6) 



*We assume systems can only exchange energy. 
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Equation (jA.6[) provides the basic relation in the canonical ensemble theory. Once T is known the other 
thermodynamic quantities can be easily obtained. 

Note that wc have chosen the average energy to be the energy of the system in this theory. One can 
show the total energy that we associate to the system on micro canonical ensemble theory corresponds to 
the most probable energy of the system, i.e. the energy level that maximizes Pr at a given temperature T . 
In practice, i.e. in the thcrmodynamical limit iV — > oo, it can be shown that these energies are equal and 
thus these two smilingly different approaches are the same. 

Finally, note that 

_ g-r cyiY>{-Er/kBT) _ exp[-(g,,. - kBTh vgr) /ksT] _ cxp{-Tr/kBT) ^ ^ 

" " Q{T, T) " Q(r, T) 

where we use S = ks Infi, which is justified by the equivalence of the two ensemble theories. Equation (jA.7[) 
shows that to maximize Pr at a fixed temperature, wc need to minimize J> over all admissible states r. To 
summarize, we have shown that minimizing the Helmholtz free energy at a temperature T (and constant 
volume) is equivalent to finding the most probable energy level, which is the total energy of the system. 
Note that this minimization should be done over all variables that determine the free energy. 



Q(T,T) 
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